function [model,errorMes] = solveQTSM(alpha,beta,psi,mu,phi,sigma,setup)

%% Allocating memory 
maxMat     = max(setup.matSelect);            % max maturity
ny         = size(setup.matSelect,2);         % number of bond yields
nxtil      = 2*setup.nm + setup.nn;           % number of pricing factors
errorMes   = 0;
% Law of motion for factors under Q 
h0Q = phi*mu;
hxQ = eye(nxtil)-phi;

% Check if co-variance matrix for innovations is positive semi-definite
sigma2 = sigma*sigma';
if min(eig(sigma2))<1D-14
    errorMes = 1;
    model = [];
    return
end
%% Compute bond pricing coefficients
% Initialise memory for bond price coefficients in state 1
A1 = zeros(1,maxMat);
B1 = zeros(nxtil,maxMat);
C1 = zeros(nxtil,nxtil,maxMat);

% Compute recursive coefficients for bond prices
for n = 1:maxMat   
    % One-period bond prices use the initial conditions
    if n == 1
        % Compute the recursive bond pricing coefficients.  All terms
        % involving multiplication by A0, B0 or C0 are zero.
        A1(1,n) = - alpha;
        B1(:,n) = - beta;
        C1(:,:,n) = - psi;

    % Longer bond prices depend on the previous maturity
    else
        gamma    = eye(nxtil)-2*sigma'*C1(:,:,n-1)*sigma;
        %invgamma = inv(gamma);
        invgamma = gamma\eye(nxtil);
        %A1(1,n) = A1(1,n-1) - alpha + B1(:,n-1)'*phi*mu ...
        %    + mu'*phi'*C1(:,:,n-1)*phi*mu ...
        %    + 0.5*B1(:,n-1)'*sigma*invgamma*sigma'*B1(:,n-1) ...
        %    + 2*(phi*mu)'*C1(:,:,n-1)*sigma*invgamma*sigma'*B1(:,n-1)...
        %    + 2*(phi*mu)'*C1(:,:,n-1)*sigma*invgamma*sigma'*C1(:,:,n-1)'*(phi*mu) ...
        %    - 0.5*log(det(gamma)); 
        detGamma = det(gamma);
        if detGamma < 0
            errorMes = 1;
            model = struct();
            return
        end
        A1(1,n) = A1(1,n-1) - alpha + B1(:,n-1)'*h0Q ...
            + h0Q'*C1(:,:,n-1)*h0Q ...
            + 0.5*B1(:,n-1)'*sigma*invgamma*sigma'*B1(:,n-1) ...
            + 2*h0Q'*C1(:,:,n-1)*sigma*invgamma*sigma'*B1(:,n-1)...
            + 2*h0Q'*C1(:,:,n-1)*sigma*invgamma*sigma'*C1(:,:,n-1)'*h0Q ...
            - 0.5*log(detGamma);       
        
        %B1(:,n) = (B1(:,n-1)'*(eye(nxtil)-phi)-beta'...
        %            +2*mu'*phi'*C1(:,:,n-1)*(eye(nxtil)-phi)...
        %            +2*B1(:,n-1)'*sigma*invgamma*sigma'*C1(:,:,n-1)'*(eye(nxtil)-phi)...
        %            +4*(phi*mu)'*C1(:,:,n-1)*sigma*invgamma*sigma'*C1(:,:,n-1)'*(eye(nxtil)-phi))';
        B1(:,n) = (B1(:,n-1)'*(eye(nxtil)-phi)-beta'...
                    +2*h0Q'*C1(:,:,n-1)*hxQ ...
                    +2*B1(:,n-1)'*sigma*invgamma*sigma'*C1(:,:,n-1)'*hxQ...
                    +4*h0Q'*C1(:,:,n-1)*sigma*invgamma*sigma'*C1(:,:,n-1)'*hxQ)';     

        %C1(:,:,n) = (eye(nxtil)-phi)'*C1(:,:,n-1)*(eye(nxtil)-phi) - psi ...
        %            + 2*(eye(nxtil)-phi)*C1(:,:,n-1)*sigma*invgamma*sigma'*C1(:,:,n-1)'*(eye(nxtil)-phi);
        C1(:,:,n) = hxQ'*C1(:,:,n-1)*hxQ - psi ...
                    + 2*hxQ*C1(:,:,n-1)*sigma*invgamma*sigma'*C1(:,:,n-1)'*hxQ;

    end
end


%% Reporting output in the desired format: annualized yields in pct.
% Spot: y_{n,t} = -(1/n)*p_{n,t}
nx         = nxtil;

% The short rate
r0   = 12*(-A1(1)/1);
rx   = 12*(-B1(:,1)'/1);
rxx  = 12*(-reshape(C1(:,:,1),1,nx^2)/1);
g0   = 12*(-A1(1,setup.matSelect)'./setup.matSelect');
gx   = 12*(-B1(:,setup.matSelect)'./repmat(setup.matSelect',1,nx));
gxx  = zeros(ny,nx^2);
for i=1:ny
    gxx(i,:) = 12*(-reshape(C1(:,:,setup.matSelect(1,i)),nx^2,1)./repmat(setup.matSelect(1,i),1,nx^2)');
end

% Storing output
model = struct('ny',length(g0),'nx',nx,'g0',g0,'gx',gx,'gxx',gxx,'h0Q',h0Q,'hxQ',hxQ,'matSelect',setup.matSelect,...
    'r0',r0,'rx',rx,'rxx',rxx,...
    'alpha',alpha,'beta',beta,'psi',psi,'phi',phi,'mu',mu,'sigmaxtil',sigma,...
    'A',A1,'B',B1,'C',C1);


end